library(tidyverse) # Data wrangling and visualization
library(colorspace) # Color manipulation
library(cmdstanr) # Interface to 'CmdStan'
library(brms) # Bayesian regression models
library(rstan) # Interface to Stan
library(ggeffects) # Marginal effects plots
library(DHARMa) # Residual diagnostics
library(emmeans) # Estimated marginal means
library(tidybayes) # Tidy data for Bayesian models
library(vegan) # Community ecology analysis
library(EcolUtils) # Ecology utilities
library(patchwork) # Combine ggplots
library(HDInterval) # HPD intervals
library(report) # Automated reporting
library(ggsignif) # Significance levels in ggplot
library(ggimage) # Images in ggplot
library(cowplot) # Plot composition in ggplot
library(paletteer) # For colour palettes
library(pilot) # For colour palettesLoad packages
Load packages needed for data manipulation and analyses.
Functions
Custom ggplot theme for visualisation.
my.theme <- function(){
theme_classic() +
theme(text = element_text(family = "Avenir Next"),
axis.title.y = element_text(margin = margin(t = 0,r = 20,b = 0,l = 0)),
axis.title.x = element_text(margin = margin(t = 20,r = 0,b = 0,l = 0)),
plot.margin = unit(c(5, 10, 5, 10), units = "mm"),
strip.background = element_rect(fill = "#CCCCFF"),
strip.text.x = element_text(size = 20),
axis.title = element_text(size = 20),
axis.text = element_text(size = 18),
legend.text = element_text(size = 15),
legend.title = element_text(size = 15))
}Load data
All mammal data
mammal_data <- read.csv("mammal_data.csv")mammal_data_vocal <- mammal_data %>%
filter(vocal == "yes",
!common.name %in% c("Yellow-bellied Sheathtail-bat")) # Removed because mostly out of human audible rangeFunction to summarise the total number of species per method for each plot across all sites. Adding missing zero values for methods that didn’t capture any species at a given plot.
summarise_mammal_data <- function(df) {
site_plots <- c("Tarcutta.DryA", "Tarcutta.DryB", "Tarcutta.WetA", "Tarcutta.WetB",
"Duval.DryA", "Duval.DryB", "Duval.WetA", "Duval.WetB",
"Mourachan.DryA", "Mourachan.DryB", "Mourachan.WetA", "Mourachan.WetB",
"Wambiana.DryA", "Wambiana.DryB", "Wambiana.WetA", "Wambiana.WetB",
"Undara.DryA", "Undara.DryB", "Undara.WetA", "Undara.WetB",
"Rinyirru.DryA", "Rinyirru.DryB", "Rinyirru.WetA", "Rinyirru.WetB")
sites <- c("Tarcutta", "Duval", "Mourachan", "Wambiana", "Undara", "Rinyirru")
methods <- c("PAM.all", "OBM", "PAM.survey", "camera")
df <- df %>%
unite(site.plot, c(site, plot), sep = ".", remove = FALSE) %>%
select(-plot) %>%
group_by(site, site.plot, assessment.method) %>%
summarise(richness = length(unique(common.name))) %>% # Needs to be common.name because of Dingo and Dog having the same scientific.name
ungroup() %>%
select(site, site.plot, assessment.method, richness) %>%
full_join(crossing(site.plot = unique(.$site.plot),
assessment.method = unique(.$assessment.method)),
by = c("site.plot", "assessment.method")) %>%
replace_na(list(richness = 0)) %>%
mutate(site.plot = factor(site.plot, levels = site_plots),
site = str_extract(site.plot, "^[^.]*"),
site = factor(site, levels = sites),
assessment.method = factor(assessment.method, levels = methods),
richness = as.numeric(richness))
return(df)
}Run function.
mammal_data_summary <- summarise_mammal_data(mammal_data)
mammal_data_vocal_summary <- summarise_mammal_data(mammal_data_vocal)PAM data
pam_data <- read_csv("mammals_PAM_daily_all.csv", show_col_types = FALSE) %>%
filter(MANUAL.ID != "UNSURE") %>%
mutate(MANUAL.ID = ifelse(MANUAL.ID == TRUE, 1, 0)) %>%
group_by(template) %>%
filter(any(MANUAL.ID == 1)) %>% #Filter for templates with at least one true positive
ungroup()BirdNET Embeddings Performance Check
Custom function
Creating a smaller data frame specifically for plotting a smooth line (100 points per template) in plots.
# Function to return predictions over extended range
run_model <- function(df) {
model_logistic <- glm(MANUAL.ID ~ eucledian.distance,
family = binomial(link = "logit"),
data = df)
# Create extended x range
new_x <- data.frame(
eucledian.distance = seq(min(df$eucledian.distance),
max(df$eucledian.distance),
length.out = 100)
)
# Get predictions and SE
preds <- predict(model_logistic, newdata = new_x, type = "link", se.fit = TRUE)
# Add predictions to new_x
new_x$fit <- plogis(preds$fit)
new_x$lower <- plogis(preds$fit - 1.96 * preds$se.fit)
new_x$upper <- plogis(preds$fit + 1.96 * preds$se.fit)
# Add template column to predictions
new_x$template <- unique(df$template)
return(new_x) # Return only the predictions dataframe
}
# Get the predictions
model_predictions <- pam_data %>%
group_by(template) %>%
group_modify(~run_model(.x)) %>%
ungroup()Recall calculations Embeddings
Sample from actual values (does not assume normality).
sample.vec <- function(x, size = 1) {
x[sample.int(length(x), size = size, replace = TRUE)]
}
# Function that uses a simple for-loop to compare random draws from positives vs negatives
calculate_probability_loop <- function(distancesPositives, distancesNegatives, nComparisons = 10000) {
set.seed(123)
# Set up a vector to record comparisons
comparisons <- c()
for (i in seq_len(nComparisons)) {
# Draw one random value from positives and one from negatives
comparisons <- c(
comparisons,
sample.vec(distancesPositives, 1) < sample.vec(distancesNegatives, 1)
)
}
# Return the fraction of times pos < neg
sum(comparisons) / length(comparisons)
}
# Now apply this by template
template_probabilities <- pam_data %>%
group_by(template) %>%
summarise(
probability = calculate_probability_loop(
eucledian.distance[MANUAL.ID == 1],
eucledian.distance[MANUAL.ID == 0],
10000
)
)Visualisation
(combined_plot <- ggplot(pam_data, aes(x = eucledian.distance, y = MANUAL.ID)) +
# Add CI ribbon
geom_ribbon(data = model_predictions,
aes(y = fit, ymin = lower, ymax = upper),
fill = "grey70",
alpha = 0.3) +
# Add CI boundary lines
geom_line(data = model_predictions,
aes(y = lower),
linetype = "dashed",
color = "black",
linewidth = 0.8) +
geom_line(data = model_predictions,
aes(y = upper),
linetype = "dashed",
color = "black",
linewidth = 0.8) +
# Add mean prediction line
geom_line(data = model_predictions,
aes(y = fit),
color = "black",
linewidth = 2) +
geom_rug(data = subset(pam_data, MANUAL.ID == 0),
aes(x = eucledian.distance),
color = "#5093CB", alpha = 0.3,
length = unit(0.05, "npc"),
linewidth = 0.8,
sides = "b",
inherit.aes = FALSE) +
geom_rug(data = subset(pam_data, MANUAL.ID == 1),
aes(x = eucledian.distance),
color = "#F78303", alpha = 0.3,
length = unit(0.05, "npc"),
linewidth = 0.8,
sides = "t",
inherit.aes = FALSE) +
geom_text(data = template_probabilities,
aes(x = Inf, y = Inf, label = sprintf("%.2f", probability)),
hjust = 1.3, vjust = 2, inherit.aes = FALSE) +
facet_wrap(~ template, scales = "free_x", ncol = 10) +
labs(x = "Euclidean Distance",
y = "Probability of correct prediction") +
scale_y_continuous(limits = c(-0.1, 1.1), breaks = seq(0, 1, 0.25)) +
theme_bw() +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.title = element_text(size = 23, face = "bold"),
axis.title.x = element_text(margin = margin(t = 25, r = 0, b = 0, l = 0)),
axis.title.y = element_text(margin = margin(t = 0, r = 25, b = 0, l = 0)),
plot.margin = unit(c(5, 10, 5, 10), units = "mm")))
Mammal Traits Effect
Imported from COMBINE (https://doi.org/10.1002/ecy.3344) and added five species not in the database using Field Companion to the Mammals of Australia and HomeRangeData (https://doi.org/10.1111/geb.13625). HomeRangeData was also used to fill in missing home range information. Detectability was calculated by divided true positive detections by overall detections and added.
trait_detect <- read_csv("trait_detect.csv") %>%
mutate(fossoriality = as.factor(fossoriality),
activity_cycle = as.factor(activity_cycle),
terrestrial_volant = as.factor(terrestrial_volant))Rows: 19 Columns: 8
── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): scientific.name
dbl (7): adult_mass_g, dispersal_km, fossoriality, home_range_km2, activity_cycle, terrestrial_volant, detectability
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.custom_labels <- function(x) {
ifelse(x >= 1,
scales::comma(x, accuracy = 1), # Whole numbers for ≥ 1
ifelse(x >= 0.1,
scales::comma(x, accuracy = 0.1), # 1 decimal for 0.1 to < 1
scales::comma(x, accuracy = 0.01))) # 2 decimals for < 0.1
}
p1 <- trait_detect %>%
ggplot(aes(adult_mass_g/1000, detectability)) +
geom_smooth(method = "lm", col = "#F78303", linewidth = 1.5) +
geom_point(size = 2.5, alpha = 0.8) +
my.theme() +
scale_y_continuous(name = "Detection probability",
breaks = seq(0, 0.10, by = 0.05)) +
coord_cartesian(ylim = c(-0.02, 0.12)) +
scale_x_log10(name = "Adult mass (kg)",
labels = custom_labels,
expand = expansion(mult = c(-0.14, 0.05)))p2 <- ggplot(trait_detect, aes(dispersal_km, detectability)) +
geom_smooth(method = "lm", col = "#F78303", linewidth = 1.5) +
geom_point(size = 2.5, alpha = 0.8) +
my.theme() +
scale_y_continuous(name = "Detection probability",
breaks = seq(0, 0.10, by = 0.05)) +
coord_cartesian(ylim = c(-0.02, 0.12)) +
scale_x_continuous(name = "Dispersal (km)",
breaks = seq(0,16,by=2),
expand = expansion(mult = c(0.02, 0.05)))p3 <- trait_detect %>%
ggplot(aes(home_range_km2, detectability)) +
geom_smooth(method = "lm", col = "#F78303", linewidth = 1.5) +
geom_point(size = 2.5, alpha = 0.8) +
my.theme() +
scale_y_continuous(name = "Detection probability",
breaks = seq(0, 0.10, by = 0.05)) +
coord_cartesian(ylim = c(-0.02, 0.12)) +
scale_x_log10(name = "Home range (km)",
labels = custom_labels,
expand = expansion(mult = c(-0.12, 0.05)))pal <- c("#204366","#CB8D79")
pal_dark <- after_scale(darken(pal, .1, space = "HLS"))
p4 <- trait_detect %>%
drop_na(fossoriality) %>%
ggplot(aes(fossoriality, detectability, fill = fossoriality)) +
geom_boxplot(aes(color = fossoriality,
fill = after_scale(desaturate(lighten(color, .8), .4))),
position = position_dodge2(width = 0.4),
width = .6, outlier.shape = NA) +
geom_point(aes(fill = fossoriality), color = "transparent",
position = position_dodge2(width = 0.4),
shape = 21, stroke = .4, size = 3.5, alpha = .4) +
my.theme() +
scale_fill_manual(values = pal, guide = "none") +
scale_colour_manual(values = pal_dark,guide = "none") +
scale_y_continuous(name = "Detection probability",
breaks = seq(0, 0.10, by = 0.05)) +
scale_x_discrete(name = "Fossoriality",
labels = c("Fossorial", "Above ground dwelling"))pal2 <- c("#204366","#CB8D79", "#324742")
pal2_dark <- after_scale(darken(pal2, .1, space = "HLS"))
p5 <- trait_detect %>%
drop_na(activity_cycle) %>%
ggplot(aes(activity_cycle, detectability, fill = activity_cycle)) +
geom_boxplot(aes(color = activity_cycle,
fill = after_scale(desaturate(lighten(color, .8), .4))),
position = position_dodge2(width = 0.4),
width = .6, outlier.shape = NA) +
geom_point(aes(fill = activity_cycle), color = "transparent",
position = position_dodge2(width = 0.4),
shape = 21, stroke = .4, size = 3.5, alpha = .4) +
my.theme() +
scale_fill_manual(values = pal2, guide = "none") +
scale_colour_manual(values = pal2_dark,guide = "none") +
scale_y_continuous(name = "Detection probability",
breaks = seq(0, 0.10, by = 0.05)) +
scale_x_discrete(name = "Activity cycle",
labels = c("Nocturnal", "Cathemeral", "Diurnal"))pal3 <- c("#204366","#CB8D79")
pal3_dark <- after_scale(darken(pal3, .1, space = "HLS"))
p6 <- trait_detect %>%
drop_na(terrestrial_volant) %>%
ggplot(aes(terrestrial_volant, detectability, fill = terrestrial_volant)) +
geom_boxplot(aes(color = terrestrial_volant,
fill = after_scale(desaturate(lighten(color, .8), .4))),
position = position_dodge2(width = 0.4),
width = .6, outlier.shape = NA) +
geom_point(aes(fill = terrestrial_volant), color = "transparent",
position = position_dodge2(width = 0.4),
shape = 21, stroke = .4, size = 3.5, alpha = .4) +
my.theme() +
scale_fill_manual(values = pal3, guide = "none") +
scale_colour_manual(values = pal3_dark,guide = "none") +
scale_y_continuous(name = "Detection probability",
breaks = seq(0, 0.10, by = 0.05)) +
scale_x_discrete(name = "Volant",
labels = c("Non-flying", "Flying"))(traits_plot <- p1 + p2 + p3 + p4 + p5 + p6 +
plot_annotation(tag_levels = "A") &
theme(plot.tag = element_text(size = 20, face = "bold")))
Investigation
trait_detect_clean <- trait_detect %>%
drop_na(detectability)
model1 <- glm(detectability ~ adult_mass_g, data = trait_detect_clean)
summary(model1)
Call:
glm(formula = detectability ~ adult_mass_g, data = trait_detect_clean)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.876e-02 8.854e-03 4.378 0.00041 ***
adult_mass_g -1.834e-08 3.461e-08 -0.530 0.60299
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.001249275)
Null deviance: 0.021589 on 18 degrees of freedom
Residual deviance: 0.021238 on 17 degrees of freedom
AIC: -69.212
Number of Fisher Scoring iterations: 2model2 <- glm(detectability ~ dispersal_km, data = trait_detect_clean)
summary(model2)
Call:
glm(formula = detectability ~ dispersal_km, data = trait_detect_clean)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.039930 0.017131 2.331 0.0481 *
dispersal_km -0.001082 0.002359 -0.459 0.6587
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.001204296)
Null deviance: 0.0098876 on 9 degrees of freedom
Residual deviance: 0.0096344 on 8 degrees of freedom
(9 observations deleted due to missingness)
AIC: -35.071
Number of Fisher Scoring iterations: 2model3 <- glm(detectability ~ fossoriality, data = trait_detect_clean)
summary(model3)
Call:
glm(formula = detectability ~ fossoriality, data = trait_detect_clean)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.039218 0.010709 3.662 0.00193 **
fossoriality2 -0.005552 0.016504 -0.336 0.74067
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.001261516)
Null deviance: 0.021589 on 18 degrees of freedom
Residual deviance: 0.021446 on 17 degrees of freedom
AIC: -69.027
Number of Fisher Scoring iterations: 2model4 <- glm(detectability ~ home_range_km2, data = trait_detect_clean)
summary(model4)
Call:
glm(formula = detectability ~ home_range_km2, data = trait_detect_clean)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.722e-02 9.998e-03 3.723 0.00256 **
home_range_km2 -1.391e-06 5.429e-06 -0.256 0.80183
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.001394961)
Null deviance: 0.018226 on 14 degrees of freedom
Residual deviance: 0.018134 on 13 degrees of freedom
(4 observations deleted due to missingness)
AIC: -52.202
Number of Fisher Scoring iterations: 2model5 <- glm(detectability ~ activity_cycle, data = trait_detect_clean)
summary(model5)
Call:
glm(formula = detectability ~ activity_cycle, data = trait_detect_clean)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.035560 0.012156 2.925 0.00991 **
activity_cycle2 0.005976 0.018378 0.325 0.74925
activity_cycle3 -0.005584 0.024312 -0.230 0.82123
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.00132988)
Null deviance: 0.021589 on 18 degrees of freedom
Residual deviance: 0.021278 on 16 degrees of freedom
AIC: -67.176
Number of Fisher Scoring iterations: 2aov1 <- aov(detectability ~ activity_cycle, data = trait_detect_clean)
summary(aov1) Df Sum Sq Mean Sq F value Pr(>F)
activity_cycle 2 0.00031 0.0001552 0.117 0.891
Residuals 16 0.02128 0.0013299 model6 <- glm(detectability ~ terrestrial_volant, data = trait_detect_clean)
summary(model6)
Call:
glm(formula = detectability ~ terrestrial_volant, data = trait_detect_clean)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.036544 0.009199 3.972 0.000984 ***
terrestrial_volant1 0.001597 0.020050 0.080 0.937460
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.001269441)
Null deviance: 0.021589 on 18 degrees of freedom
Residual deviance: 0.021580 on 17 degrees of freedom
AIC: -68.908
Number of Fisher Scoring iterations: 2Species Richness
Given the inherent inability of passive acoustic monitoring to detect non-vocalising species, the mammal community was considered in two ways: i) entire mammal community (All Mammals), and ii) vocal mammals only (Vocal Mammals).
All Mammals
Considering the entire mammal community with vocal and non-vocal species.
Bayesian Analysis
Priors
Defining priors for Bayesian models.
mammal_data_summary %>% summarise(median(log(richness)))# A tibble: 1 × 1
`median(log(richness))`
<dbl>
1 1.79mammal_data_summary %>% summarise(mad(log(richness)))# A tibble: 1 × 1
`mad(log(richness))`
<dbl>
1 0.601Priors for the intercept are 1.8 and 1.
log(sd(mammal_data_summary$richness))/apply(model.matrix(~assessment.method, data = mammal_data_summary), 2, sd) (Intercept) assessment.methodOBM assessment.methodPAM.survey assessment.methodcamera
Inf 2.829888 2.829888 2.829888 Priors for the slope are 2.8.
Fit model
priors1 <- prior(normal(1.8,1), class = "Intercept") +
prior(normal(0,2.8), class = "b") +
prior(student_t(3,0,2.8), class = "sd")
methods.form1 <- bf(richness ~ assessment.method + (1|site/site.plot),
family = poisson(link = "log"))
methods.brm1 <- brm(methods.form1,
data = mammal_data_summary,
prior = priors1,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 3,
warmup = 1000,
backend = 'cmdstanr',
save_pars = save_pars(all = TRUE))Start samplingRunning MCMC with 3 parallel chains...
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.2 seconds.Prior vs Posterior
methods.brm1b <- methods.brm1 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 3,
cores = 3,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE))The desired updates require recompiling the modelStart samplingRunning MCMC with 3 parallel chains...
Chain 2 finished in 0.9 seconds.
Chain 3 finished in 0.9 seconds.
Chain 1 finished in 1.1 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.9 seconds.
Total execution time: 1.1 seconds.Warning: 4 of 2400 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.priors2 <- prior(normal(1.8,1), class = "Intercept") +
prior(normal(0,2.8), class = "b") +
prior(student_t(3,0,2.8), class = "sd")
methods.form2 <- bf(richness ~ assessment.method + (1|site/site.plot),
family="negbinomial")
methods.brm2 <- brm(methods.form2,
data = mammal_data_summary,
prior = priors2,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 8,
warmup = 1000,
backend = 'cmdstanr',
save_pars = save_pars(all = TRUE))Start samplingRunning MCMC with 3 parallel chains...
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.2 seconds.
Chain 3 finished in 0.2 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.3 seconds.Warning: 1 of 2400 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.Prior vs Posterior
methods.brm2b <- methods.brm2 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 8,
cores = 3,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE))The desired updates require recompiling the modelStart samplingRunning MCMC with 3 parallel chains...
Chain 2 finished in 1.1 seconds.
Chain 3 finished in 1.2 seconds.
Chain 1 finished in 26.8 seconds.
All 3 chains finished successfully.
Mean chain execution time: 9.7 seconds.
Total execution time: 26.9 seconds.Warning: 166 of 2400 (7.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.Warning: 712 of 2400 (30.0%) transitions hit the maximum treedepth limit of 10.
See https://mc-stan.org/misc/warnings for details.Warning: 1 of 3 chains have a NaN E-BFMI.
See https://mc-stan.org/misc/warnings for details.priors3 <- prior(normal(1.8,1), class = "Intercept") +
prior(normal(0,2.8), class = "b") +
prior(student_t(3,0,2.8), class = "sd")
methods.form3 <- bf(richness ~ assessment.method + (1|site/site.plot),
family="negbinomial2")
methods.brm3 <- brm(methods.form3,
data = mammal_data_summary,
prior = priors3,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 1,
warmup = 1000,
control = list(adapt_delta=0.99),
backend = 'cmdstanr',
save_pars = save_pars(all = TRUE))Start samplingRunning MCMC with 3 parallel chains...
Chain 1 finished in 0.3 seconds.
Chain 2 finished in 0.3 seconds.
Chain 3 finished in 0.3 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.5 seconds.Prior vs Posterior
methods.brm3b <- methods.brm3 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 2,
cores = 3,
control = list(adapt_delta=0.99),
backend = "cmdstanr",
save_pars = save_pars(all = TRUE))The desired updates require recompiling the modelStart samplingRunning MCMC with 3 parallel chains...
Chain 1 finished in 2.4 seconds.
Chain 2 finished in 2.4 seconds.
Chain 3 finished in 31.2 seconds.
All 3 chains finished successfully.
Mean chain execution time: 12.0 seconds.
Total execution time: 31.3 seconds.Warning: 800 of 2400 (33.0%) transitions hit the maximum treedepth limit of 10.
See https://mc-stan.org/misc/warnings for details.Warning: 1 of 3 chains have a NaN E-BFMI.
See https://mc-stan.org/misc/warnings for details.Compare Models
(l.1b <- methods.brm1b %>% loo())
Computed from 2400 by 96 log-likelihood matrix.
Estimate SE
elpd_loo -204.3 3.6
p_loo 5.5 0.6
looic 408.6 7.2
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.(l.2b <- methods.brm2b %>% loo())
Computed from 2400 by 96 log-likelihood matrix.
Estimate SE
elpd_loo -208.5 4.2
p_loo 6.4 0.8
looic 417.0 8.3
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.9]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 60 62.5% 3
(0.7, 1] (bad) 1 1.0% <NA>
(1, Inf) (very bad) 35 36.5% <NA>
See help('pareto-k-diagnostic') for details.(l.3b <- methods.brm3b %>% loo())
Computed from 2400 by 96 log-likelihood matrix.
Estimate SE
elpd_loo -2358.6 149.4
p_loo 2121.2 147.4
looic 4717.2 298.9
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.9]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 1 1.0% 1681
(0.7, 1] (bad) 0 0.0% <NA>
(1, Inf) (very bad) 95 99.0% <NA>
See help('pareto-k-diagnostic') for details.loo_compare(loo(methods.brm1b), loo(methods.brm2b), loo(methods.brm3b)) elpd_diff se_diff
methods.brm1b 0.0 0.0
methods.brm2b -4.2 1.5
methods.brm3b -2154.3 147.6Model methods.brm1b was selected as best model based on loo estimates.
Diagnostics
methods.brm1b$fit %>% stan_trace()'pars' not specified. Showing first 10 parameters by default.
methods.brm1b$fit %>% stan_ac()'pars' not specified. Showing first 10 parameters by default.
methods.brm1b$fit %>% stan_rhat()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
methods.brm1b$fit %>% stan_ess()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
methods.brm1b %>% pp_check(type = "dens_overlay", ndraws = 200)
set.seed(6)
preds <- posterior_predict(methods.brm1b, ndraws=250, summary=FALSE)
method.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = mammal_data_summary$richness,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE)
method.resids %>% plot()
method.resids %>% testDispersion()
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.37045, p-value < 2.2e-16
alternative hypothesis: two.sidedmethod.resids %>% testZeroInflation()
DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 0, p-value = 0.352
alternative hypothesis: two.sidedInvestigation
Methods pairwise comparison across sites
(diff.methods.avg <- methods.brm1b %>%
emmeans(~assessment.method) %>%
regrid(transform = "none") %>%
pairs() %>%
gather_emmeans_draws() %>%
mutate(f.change = exp(.value)) %>%
summarise("Average fractional change" = median(f.change),
"Lower HDI" = HDInterval::hdi(f.change)[1],
"Upper HDI" = HDInterval::hdi(f.change)[2],
"Probability of difference" = sum(.value > 0)/n())) #to see if there is any change# A tibble: 6 × 5
contrast `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
<chr> <dbl> <dbl> <dbl> <dbl>
1 OBM - PAM.survey 2.46 1.85 3.12 1
2 OBM - camera 1.44 1.16 1.77 0.999
3 PAM.all - OBM 0.863 0.691 1.04 0.0729
4 PAM.all - PAM.survey 2.13 1.62 2.72 1
5 PAM.all - camera 1.25 0.992 1.52 0.977
6 PAM.survey - camera 0.587 0.430 0.744 0 (diff.methods.avg.rev <- methods.brm1b %>%
emmeans(~assessment.method) %>%
regrid(transform = "none") %>%
pairs(reverse = TRUE) %>%
gather_emmeans_draws() %>%
mutate(f.change = exp(.value)) %>%
summarise("Average fractional change" = median(f.change),
"Lower HDI" = HDInterval::hdi(f.change)[1],
"Upper HDI" = HDInterval::hdi(f.change)[2],
"Probability of difference" = sum(.value > 0)/n())) #to see if there is any change# A tibble: 6 × 5
contrast `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
<chr> <dbl> <dbl> <dbl> <dbl>
1 OBM - PAM.all 1.16 0.932 1.39 0.927
2 PAM.survey - OBM 0.406 0.304 0.513 0
3 PAM.survey - PAM.all 0.469 0.359 0.601 0
4 camera - OBM 0.693 0.555 0.848 0.000833
5 camera - PAM.all 0.801 0.651 0.999 0.0229
6 camera - PAM.survey 1.70 1.26 2.17 1 Visualisation
pal <- c("#189F9F","#FF8C00", "#18709F", "#A034F0")
pal_dark <- after_scale(darken(pal, .1, space = "HLS"))
(methods.boxplot <- ggplot(mammal_data_summary, aes(assessment.method, richness)) +
geom_boxplot(aes(color = assessment.method,
fill = after_scale(desaturate(lighten(color, .8), .4))),
position = position_dodge2(width = 0.6),
width = .6, outlier.shape = NA) +
geom_point(aes(fill = assessment.method), color = "transparent",
position = position_dodge2(width = 0.6),
shape = 21, stroke = .4, size = 3.5, alpha = .3) +
my.theme() +
scale_y_continuous(limits = c(0,20),breaks = seq(0, 20, by = 2),
name = "Total Species Richness") +
scale_x_discrete(name = "", labels = c("", "", "", "")) +
scale_fill_manual(values = pal, guide = "none") +
scale_colour_manual(values = pal_dark, name = "", labels = c("PAM (all audio)","OBM", "PAM (survey period)","Camera Trap")) +
theme(legend.text = element_text(size = 14, margin = margin(t = 5)),
legend.position = c(y=0.45, x=-0.15),
legend.background = element_rect(fill = "transparent"),
plot.margin = unit(c(5, 10, 20, 10), units = "mm")) +
guides(colour = guide_legend(nrow = 1)) +
geom_signif(comparisons = list(c("PAM.all", "PAM.survey"),
c("OBM", "PAM.survey"),
c("OBM", "camera"),
c("PAM.survey", "camera")),
annotations = c("2.13, HDI 95%[1.62,2.72]",
"2.46, HDI 95%[1.85,3.12]",
"1.44, HDI 95%[1.16,1.77]",
"0.59, HDI 95%[0.43,0.74]"),
y_position = c(16, 17.5, 19, 12)))
Vocal Mammals
Considering vocal mammals only and excluding non-vocal mammals from the analysis.
Bayesian Analysis
Priors
Defining priors for Bayesian models.
mammal_data_vocal_summary %>% summarise(median(log(richness)))# A tibble: 1 × 1
`median(log(richness))`
<dbl>
1 1.61mammal_data_vocal_summary %>% summarise(mad(log(richness)))# A tibble: 1 × 1
`mad(log(richness))`
<dbl>
1 0.598Priors for the intercept are 1.6 and 1.
log(sd(mammal_data_vocal_summary$richness))/apply(model.matrix(~assessment.method, data = mammal_data_vocal_summary), 2, sd) (Intercept) assessment.methodOBM assessment.methodPAM.survey assessment.methodcamera
Inf 2.220652 2.220652 2.220652 Priors for the slope are 2.5.
Fit model
priors_vocal1 <- prior(normal(1.6,1), class = "Intercept") +
prior(normal(0,2.5), class = "b") +
prior(student_t(3,0,2.5), class = "sd")
methods.vocal.form1 <- bf(richness ~ assessment.method + (1|site/site.plot),
family = poisson(link = "log"))
methods.vocal.brm1 <- brm(methods.vocal.form1,
data = mammal_data_vocal_summary,
prior = priors_vocal1,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 3,
warmup = 1000,
backend = 'cmdstanr',
save_pars = save_pars(all = TRUE))Start samplingRunning MCMC with 3 parallel chains...
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.2 seconds.Prior vs Posterior
methods.vocal.brm1b <- methods.vocal.brm1 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 3,
cores = 3,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE))The desired updates require recompiling the modelStart samplingRunning MCMC with 3 parallel chains...
Chain 1 finished in 0.9 seconds.
Chain 2 finished in 0.9 seconds.
Chain 3 finished in 0.9 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.9 seconds.
Total execution time: 1.0 seconds.Warning: 1 of 2400 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.priors_vocal2 <- prior(normal(1.6,1), class = "Intercept") +
prior(normal(0,2.5), class = "b") +
prior(student_t(3,0,2.5), class = "sd")
methods.vocal.form2 <- bf(richness ~ assessment.method + (1|site/site.plot),
family="negbinomial")
methods.vocal.brm2 <- brm(methods.vocal.form2,
data = mammal_data_vocal_summary,
prior = priors_vocal2,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 2,
warmup = 1000,
backend = 'cmdstanr',
save_pars = save_pars(all = TRUE))Start samplingRunning MCMC with 3 parallel chains...
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.3 seconds.Warning: 2 of 2400 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.Prior vs Posterior
methods.vocal.brm2b <- methods.vocal.brm2 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 2,
cores = 3,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE))The desired updates require recompiling the modelStart samplingRunning MCMC with 3 parallel chains...
Chain 1 finished in 24.4 seconds.
Chain 2 finished in 26.0 seconds.
Chain 3 finished in 27.6 seconds.
All 3 chains finished successfully.
Mean chain execution time: 26.0 seconds.
Total execution time: 27.7 seconds.Warning: 477 of 2400 (20.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.Warning: 1923 of 2400 (80.0%) transitions hit the maximum treedepth limit of 10.
See https://mc-stan.org/misc/warnings for details.Warning: 3 of 3 chains have a NaN E-BFMI.
See https://mc-stan.org/misc/warnings for details.priors_vocal3 <- prior(normal(1.6,1), class = "Intercept") +
prior(normal(0,2.5), class = "b") +
prior(student_t(3,0,2.5), class = "sd")
methods.vocal.form3 <- bf(richness ~ assessment.method + (1|site/site.plot),
family="negbinomial2")
methods.vocal.brm3 <- brm(methods.vocal.form3,
data = mammal_data_vocal_summary,
prior = priors_vocal3,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 3,
warmup = 1000,
control = list(adapt_delta=0.99),
backend = 'cmdstanr',
save_pars = save_pars(all = TRUE))Start samplingRunning MCMC with 3 parallel chains...
Chain 1 finished in 0.3 seconds.
Chain 3 finished in 0.3 seconds.
Chain 2 finished in 0.3 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.5 seconds.Prior vs Posterior
methods.vocal.brm3b <- methods.vocal.brm3 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 3,
cores = 3,
control = list(adapt_delta=0.99),
backend = "cmdstanr",
save_pars = save_pars(all = TRUE))The desired updates require recompiling the modelStart samplingRunning MCMC with 3 parallel chains...
Chain 2 finished in 1.3 seconds.
Chain 3 finished in 1.6 seconds.
Chain 1 finished in 2.5 seconds.
All 3 chains finished successfully.
Mean chain execution time: 1.8 seconds.
Total execution time: 2.6 seconds.Warning: 715 of 2400 (30.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.Warning: 2 of 3 chains had an E-BFMI less than 0.3.
See https://mc-stan.org/misc/warnings for details.Compare Models
(l.1b <- methods.vocal.brm1b %>% loo())
Computed from 2400 by 96 log-likelihood matrix.
Estimate SE
elpd_loo -185.0 3.5
p_loo 4.3 0.5
looic 370.0 7.0
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.(l.2b <- methods.vocal.brm2b %>% loo())
Computed from 2400 by 96 log-likelihood matrix.
Estimate SE
elpd_loo -188434938.1 106228234.4
p_loo 188434718.3 106228235.0
looic 376869876.1 212456468.7
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.0]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 0 0.0% <NA>
(0.7, 1] (bad) 0 0.0% <NA>
(1, Inf) (very bad) 96 100.0% <NA>
See help('pareto-k-diagnostic') for details.(l.3b <- methods.vocal.brm3b %>% loo())
Computed from 2400 by 96 log-likelihood matrix.
Estimate SE
elpd_loo -184.6 3.4
p_loo 3.6 0.4
looic 369.2 6.8
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.9]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 79 82.3% 5
(0.7, 1] (bad) 2 2.1% <NA>
(1, Inf) (very bad) 15 15.6% <NA>
See help('pareto-k-diagnostic') for details.loo_compare(loo(methods.vocal.brm1b), loo(methods.vocal.brm2b), loo(methods.vocal.brm3b)) elpd_diff se_diff
methods.vocal.brm3b 0.0 0.0
methods.vocal.brm1b -0.4 0.4
methods.vocal.brm2b -188434753.5 106228235.0Model methods.vocal.brm1b was selected as best model based on loo estimates.
Diagnostics
methods.vocal.brm1b$fit %>% stan_trace()'pars' not specified. Showing first 10 parameters by default.
methods.vocal.brm1b$fit %>% stan_ac()'pars' not specified. Showing first 10 parameters by default.
methods.vocal.brm1b$fit %>% stan_rhat()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
methods.vocal.brm1b$fit %>% stan_ess()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
methods.vocal.brm1b %>% pp_check(type = "dens_overlay", ndraws = 200)
set.seed(6)
preds.vocal <- posterior_predict(methods.vocal.brm1b, ndraws=250, summary=FALSE)
method.resids <- createDHARMa(simulatedResponse = t(preds.vocal),
observedResponse = mammal_data_vocal_summary$richness,
fittedPredictedResponse = apply(preds.vocal, 2, median),
integerResponse = TRUE)
method.resids %>% plot()
method.resids %>% testDispersion()
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.31072, p-value < 2.2e-16
alternative hypothesis: two.sidedmethod.resids %>% testZeroInflation()
DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 0.62735, p-value = 0.784
alternative hypothesis: two.sidedInvestigation
Methods pairwise comparison across sites
(diff.methods.vocal.avg <- methods.vocal.brm1b %>%
emmeans(~assessment.method) %>%
regrid(transform = "none") %>%
pairs() %>%
gather_emmeans_draws() %>%
mutate(f.change = exp(.value)) %>%
summarise("Average fractional change" = median(f.change),
"Lower HDI" = HDInterval::hdi(f.change)[1],
"Upper HDI" = HDInterval::hdi(f.change)[2],
"Probability of difference" = sum(.value > 0)/n())) #to see if there is any change# A tibble: 6 × 5
contrast `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
<chr> <dbl> <dbl> <dbl> <dbl>
1 OBM - PAM.survey 1.34 0.989 1.74 0.980
2 OBM - camera 1.46 1.05 1.87 0.995
3 PAM.all - OBM 1.58 1.24 1.98 1
4 PAM.all - PAM.survey 2.12 1.61 2.69 1
5 PAM.all - camera 2.29 1.73 2.93 1
6 PAM.survey - camera 1.08 0.785 1.43 0.697(diff.methods.vocal.avg.rev <- methods.vocal.brm1b %>%
emmeans(~assessment.method) %>%
regrid(transform = "none") %>%
pairs(reverse = TRUE) %>%
gather_emmeans_draws() %>%
mutate(f.change = exp(.value)) %>%
summarise("Average fractional change" = median(f.change),
"Lower HDI" = HDInterval::hdi(f.change)[1],
"Upper HDI" = HDInterval::hdi(f.change)[2],
"Probability of difference" = sum(.value > 0)/n())) #to see if there is any change# A tibble: 6 × 5
contrast `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
<chr> <dbl> <dbl> <dbl> <dbl>
1 OBM - PAM.all 0.633 0.498 0.799 0
2 PAM.survey - OBM 0.746 0.555 0.974 0.0204
3 PAM.survey - PAM.all 0.472 0.357 0.600 0
4 camera - OBM 0.686 0.514 0.905 0.005
5 camera - PAM.all 0.436 0.325 0.550 0
6 camera - PAM.survey 0.923 0.681 1.23 0.303 Visualisation
pal <- c("#189F9F","#FF8C00", "#18709F","#A034F0")
pal_dark <- after_scale(darken(pal, .1, space = "HLS"))
(methods.boxplot.vocal <- ggplot(mammal_data_vocal_summary, aes(assessment.method, richness)) +
geom_boxplot(aes(color = assessment.method,
fill = after_scale(desaturate(lighten(color, .8), .4))),
position = position_dodge2(width = 0.6),
width = .6, outlier.shape = NA) +
geom_point(aes(fill = assessment.method), color = "transparent",
position = position_dodge2(width = 0.6),
shape = 21, stroke = .4, size = 3.5, alpha = .3) +
my.theme() +
scale_y_continuous(limits = c(0,20),breaks = seq(0, 20, by = 2),
name = "") +
scale_x_discrete(name = "", labels = c("", "", "", "")) +
scale_fill_manual(values = pal, guide = "none") +
scale_colour_manual(values = pal_dark, name = "", labels = c("PAM (all audio)","OBM", "PAM (survey period)","Camera Trap")) +
theme(legend.text = element_text(size = 14, margin = margin(t = 5)),
legend.position = c(y=0.45, x=-0.15),
legend.background = element_rect(fill = "transparent"),
plot.margin = unit(c(5, 10, 20, 10), units = "mm")) +
guides(colour = guide_legend(nrow = 1)) +
geom_signif(comparisons = list(c("PAM.all", "OBM"),
c("PAM.all", "PAM.survey"),
c("PAM.all", "camera"),
c("OBM", "camera")),
annotations = c("1.58, HDI 95%[1.24,1.98]",
"2.12, HDI 95%[1.61,2.69]",
"2.29, HDI 95%[1.73,2.93]",
"1.46, HDI 95%[1.05,1.87]"),
y_position = c(12.5, 14, 15.5, 9)))
Combined Figure
(methods.vocal.combined <- (methods.boxplot + ggtitle('All mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(methods.boxplot.vocal + ggtitle('Vocal mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
theme(legend.position = c(-0.25,-0.2), legend.text = element_text(size = 22),
plot.margin = unit(c(5, 10, 30, 10), units = "mm")))
Correlation PAM vocal mammals & OBM mammals
Investigate a correlation between species richness of vocal mammals detected by PAM and species richness of the entire mammal community detected by OBM.
cor_mammal_data <- mammal_data_summary |>
filter(assessment.method %in% c("OBM", "PAM.all")) |>
group_by(site.plot) |>
pivot_wider(names_from = "assessment.method", #transformed dataset from long to wide format
values_from = "richness",
values_fill = list(number=0))
(corr_plot <- ggplot(cor_mammal_data, aes(PAM.all,OBM)) +
geom_point(aes(col = site), size = 3) +
geom_smooth(method = "lm") +
my.theme() +
theme(legend.position = "bottom") +
scale_y_continuous(name = "Species Richness OBM", breaks = seq(4, 16, by = 4)) +
scale_x_continuous(name = "Species Richness PAM", breaks = seq(2, 13, by = 2)) +
scale_color_manual(name = "", values = c("seagreen", "steelblue","salmon2", "darkorchid1", "aquamarine3", "goldenrod1")) +
guides(color = guide_legend(nrow = 1)) +
annotate("text", x = 5, y = 15,
label = "R^2 == 0.21", parse = TRUE, size = 5))`geom_smooth()` using formula = 'y ~ x'
# Perform Pearson's correlation test
correlation_test <- cor.test(cor_mammal_data$PAM.all, cor_mammal_data$OBM, method = "pearson")
print(correlation_test)
Pearson's product-moment correlation
data: cor_mammal_data$PAM.all and cor_mammal_data$OBM
t = 2.4133, df = 22, p-value = 0.02458
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.06636028 0.72677443
sample estimates:
cor
0.4575097 (R.squared <- 0.4575097^2)[1] 0.2093151Community Composition
All Mammals
Considering the entire mammal community with vocal and non-vocal species.
Prepare data
mammal_community3 <- mammal_data %>%
filter(!recapture %in% c("y")) %>% # Removed recaptures
select(site, assessment.method,common.name) %>% # Selected relevant grouping variables
group_by(across(everything())) %>% # Grouped by all available columns
mutate(number = n()) %>% # Created a numbers column that counts the number of duplicate observations
ungroup() %>%
unique() %>% # Merged the duplicate rows
pivot_wider(names_from = "common.name", # Transformed dataset from long to wide format
values_from = "number",
values_fill = list(number=0)) %>% # Replaced NAs with 0
arrange(site) %>% # Ordered by site
mutate(assessment.method = factor(assessment.method)) # Turned asssessment methods into factorsIn preparation for community analyses, only keep columns and rows with values above 0.
all <- mammal_community3 %>% select(where(~ any(. != 0))) %>% filter(rowSums(.[,3:ncol(.)]) > 0)Transformation of the data to binary (presence/absence) for Jaccard dissimilarity.
all_jacc <- all %>%
select(where(is.numeric)) %>% # Removed columns with non-numerical values
decostand(method="pa", MARGIN = 1) %>% # Standardised presence/absence (method = "pa") by rows (MARGIN = 1)
cbind(all[,1:2]) %>%
select(site, assessment.method, everything())NMDS Ordination
Jaccard dissimilarity
set.seed(11)
all_nmds_jacc <- metaMDS(all_jacc[4:ncol(all_jacc)], distance="jaccard", k=3,trymax=100, autotransform = FALSE)Run 0 stress 0.1055242
Run 1 stress 0.1207629
Run 2 stress 0.1047353
... New best solution
... Procrustes: rmse 0.02496561 max resid 0.0873796
Run 3 stress 0.105525
Run 4 stress 0.1083786
Run 5 stress 0.109782
Run 6 stress 0.104735
... New best solution
... Procrustes: rmse 0.0002099691 max resid 0.0007380798
... Similar to previous best
Run 7 stress 0.1047354
... Procrustes: rmse 0.0002734112 max resid 0.0009616091
... Similar to previous best
Run 8 stress 0.114386
Run 9 stress 0.1207632
Run 10 stress 0.1133275
Run 11 stress 0.1097819
Run 12 stress 0.1123875
Run 13 stress 0.1133275
Run 14 stress 0.1145607
Run 15 stress 0.1133275
Run 16 stress 0.1055252
Run 17 stress 0.1166412
Run 18 stress 0.1145605
Run 19 stress 0.1155145
Run 20 stress 0.1200923
*** Best solution repeated 2 timesVisualisation
Jaccard dissimilarity
Function to create nMDS plots with centroids.
create_plot <- function(data, nmds_data) {
methods <- nmds_data$points
species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
df <- cbind.data.frame(data[,1:2], methods)
gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
mutate(assessment.method = factor(assessment.method, levels = c("PAM.all", "OBM","PAM.survey", "camera")))
ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
geom_point(size=3) + # Smaller outside points
geom_point(aes(x=mean.x,y=mean.y),size=5) + # This controls the centroids
scale_colour_manual(values = c("#189F9F","#FF8C00", "#18709F","#A034F0"),
name = "",
labels = c("PAM (all audio)", "OBM", "PAM (survey period)","Camera Trap")) +
geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) + # This controls the connections between points and centroids
my.theme() +
theme(legend.position = "bottom")
}
(nmds.all.plot.jacc.centroid <- create_plot(all_jacc, all_nmds_jacc))
Investigation
Adonis
Jaccard dissimilarity
all.dist2 <- vegdist(all_jacc[3:ncol(all_jacc)], 'jaccard')
(all.adonis <- adonis2(all.dist2 ~ assessment.method, data=all_jacc))Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 999
adonis2(formula = all.dist2 ~ assessment.method, data = all_jacc)
Df SumOfSqs R2 F Pr(>F)
Model 3 1.9870 0.30827 2.971 0.001 ***
Residual 20 4.4586 0.69173
Total 23 6.4456 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Permutation test: significant difference between assessment methods for mammal communities.
#Let's which method is different from which
(adonis.methods.jacc <- adonis.pair(all.dist2, all_jacc$assessment.method, corr.method = "holm")) combination SumsOfSqs MeanSqs F.Model R2 P.value P.value.corrected
1 camera <-> OBM 0.4559165 0.4559165 1.692142 0.1447247 0.105894106 0.211788212
2 camera <-> PAM.all 0.9499132 0.9499132 4.958432 0.3314807 0.000999001 0.005994006
3 camera <-> PAM.survey 1.1114297 1.1114297 4.864324 0.3272483 0.003996004 0.015984016
4 OBM <-> PAM.all 0.5417048 0.5417048 2.492023 0.1994891 0.018981019 0.056943057
5 OBM <-> PAM.survey 0.6990302 0.6990302 2.748990 0.2156241 0.002997003 0.014985015
6 PAM.all <-> PAM.survey 0.2159058 0.2159058 1.223748 0.1090320 0.290709291 0.290709291#comparing all methods with each otherVocal Mammals
Considering vocal mammals only and excluding non-vocal mammals from the analysis.
Prepare data
mammal_community4 <- mammal_data_vocal %>%
filter(!recapture %in% c("y")) %>%
select(site, assessment.method,common.name) %>% #selected relevant grouping variables
group_by(across(everything())) %>% #groups by all available columns
mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
ungroup() %>%
unique() %>% #merged the duplicate rows
pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
values_from = "number",
values_fill = list(number=0)) %>% #replaced NAs with 0
arrange(site) %>%
mutate(assessment.method = factor(assessment.method))all_vocal <- mammal_community4 %>% select(where(~ any(. != 0))) %>% filter(rowSums(.[,3:ncol(.)]) > 0)#Function decostand() standardises presence/absence (method = "pa") by rows (MARGIN = 1)
all_vocal_jacc <- all_vocal %>%
select(where(is.numeric)) %>%
decostand(method="pa", MARGIN = 1) %>%
cbind(all_vocal[,1:2]) %>%
select(site, assessment.method, everything())NMDS Ordination
Jaccard dissimilarity
set.seed(11)
all_nmds_jacc_vocal <- metaMDS(all_vocal_jacc[4:ncol(all_vocal_jacc)], distance="jaccard", k=3,trymax=100, autotransform = FALSE)Visualisation
Jaccard dissimilarity
create_plot <- function(data, nmds_data) {
methods <- nmds_data$points
species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
df <- cbind.data.frame(data[,1:2], methods)
gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
mutate(assessment.method = factor(assessment.method, levels = c("PAM.all", "OBM","PAM.survey", "camera")))
ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
geom_point(size=3) +
geom_point(aes(x=mean.x,y=mean.y),size=5) +
scale_colour_manual(values = c("#189F9F","#FF8C00", "#18709F","#A034F0"),
name = "",
labels = c("PAM (all audio)", "OBM", "PAM (survey period)","Camera Trap")) +
geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) +
my.theme() + theme(legend.position = "bottom")
}
(nmds.all.vocal.plot.jacc.centroid <- create_plot(all_vocal_jacc, all_nmds_jacc_vocal))
Investigation
Adonis
Jaccard dissimilarity
all.dist_vocal <- vegdist(all_vocal_jacc[3:ncol(all_vocal_jacc)], 'jaccard')
(all.adonis_vocal <- adonis2(all.dist_vocal ~ assessment.method, data=all_vocal_jacc))Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 999
adonis2(formula = all.dist_vocal ~ assessment.method, data = all_vocal_jacc)
Df SumOfSqs R2 F Pr(>F)
Model 3 1.3414 0.2559 2.2927 0.01 **
Residual 20 3.9006 0.7441
Total 23 5.2420 1.0000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Permutation test: significant difference between assessment methods for mammal communities.
#Let's which method is different from which
(adonis.methods.jacc_vocal <- adonis.pair(all.dist_vocal, all_vocal_jacc$assessment.method, corr.method = "holm")) combination SumsOfSqs MeanSqs F.Model R2 P.value P.value.corrected
1 camera <-> OBM 0.4621827 0.4621827 2.1635228 0.17786976 0.071928072 0.28771229
2 camera <-> PAM.all 0.6947730 0.6947730 4.1757138 0.29456815 0.004995005 0.02497502
3 camera <-> PAM.survey 0.8980116 0.8980116 4.4172857 0.30638816 0.003996004 0.02397602
4 OBM <-> PAM.all 0.1200808 0.1200808 0.6429681 0.06041248 0.713286713 0.71328671
5 OBM <-> PAM.survey 0.2919053 0.2919053 1.3050666 0.11544086 0.233766234 0.70129870
6 PAM.all <-> PAM.survey 0.2159058 0.2159058 1.2237477 0.10903200 0.290709291 0.70129870#comparing all methods with each otherCombined Figure
(nmds.jacc.combined <- (nmds.all.plot.jacc.centroid + ggtitle('All mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(nmds.all.vocal.plot.jacc.centroid + ggtitle('Vocal mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
theme(legend.position = c(-0.25,-0.2), legend.text = element_text(size = 22),
plot.margin = unit(c(5, 10, 30, 10), units = "mm")) +
guides(colour = guide_legend(nrow = 1)))
Survey Effort (Species accumulation curves)
All Mammals
Considering the entire mammal community with vocal and non-vocal species.
Prepare data
Merged duplicate rows for duplicate observations and transformed the dataset to a wide format.
mammal_community <- mammal_data %>%
select(site, SamplingDay, assessment.method,common.name) %>% #selected relevant grouping variables
group_by(across(everything())) %>% #groups by all available columns
mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
ungroup() %>%
unique() %>% #merged the duplicate rows
pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
values_from = "number",
values_fill = list(number=0)) %>% #replaced NAs with 0
arrange(site) %>%
filter(assessment.method != "PAM.all") %>%
mutate(SamplingDay = as.factor(SamplingDay))Added a category for the total richness.
result.i <- vector("list",length(unique(mammal_community$site))) #created an empty list that was filled as the loop ran
for(i in 1:length(unique(mammal_community$site))){ #looped over each study site
site.i <- mammal_community[mammal_community$site == unique(mammal_community$site)[i],] #subset to the site
result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] #subset to the SamplingDay
meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("total"))
SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
colnames(total.a) <- colnames(mammal_community)
result.a[[a]] <- total.a} #added method results to list
result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop
result.i <- do.call("rbind.data.frame",result.i)
mammal_community2 <- rbind.data.frame(mammal_community, result.i)Added a remote sensing category to assessment methods. Remote sensing was a combination of PAM and camera trapping.
result.j <- vector("list",length(unique(mammal_community$site))) #created an empty list that was filled as the loop ran
for(i in 1:length(unique(mammal_community$site))){ #looped over each study site
site.i <- mammal_community[mammal_community$site == unique(mammal_community$site)[i],] #subset to the site
result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] %>%
filter(assessment.method %in% c("camera", "PAM.survey"))
meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("remote"))
SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
colnames(total.a) <- colnames(mammal_community)
result.a[[a]] <- total.a} #added method results to list
result.j[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop
result.j <- do.call("rbind.data.frame",result.j)
mammal_community2 <- rbind.data.frame(mammal_community2, result.j) %>%
drop_na()Loop
Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.
# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 30, ncol = (ncol(mammal_community2)))) %>%
replace(is.na(.), 0) %>%
set_names(names(mammal_community2)) %>%
select(-(1:3))
result.i <- vector("list",length(unique(mammal_community2$site))) #created an empty list that was filled as the loop ran
set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(mammal_community2$site))){ #looped over each study site
site.i <- mammal_community2[mammal_community2$site == unique(mammal_community2$site)[i],] #subset to the site
dates.i <- ifelse(unique(site.i$site) %in% c("Tarcutta", "Undara", "Wambiana", "Rinyirru"), 28, 21) #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day
result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
sampling.a <- method.a$assessment.method[1] #saved the name of the method
method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
accum.a <- specaccum(method.a,"random") #calculated data
#created data frame of metadata, richness, and standard deviation
accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
result.a[[a]] <- accum.a}} #added method results to list
result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loopTransformed resulting list to a data frame.
mammals_wide.result <- do.call("rbind.data.frame",result.i) %>% #compressed result.i list into data frame
mutate(assessment.method = factor(assessment.method, levels = c("OBM",
"camera",
"PAM.survey",
"remote",
"total")),
site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))Visualisation
Plot of rarefaction curves for each method split by sites.
(specaccum_mammals <- ggplot(mammals_wide.result, aes(x = day, y = richness, col = assessment.method)) +
geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.3))),
linetype = 0.1) +
geom_vline(xintercept=c(7,14,21, 28), linetype="dotted", colour = "black") +
geom_line(aes(group = assessment.method),linewidth = 1.3) +
my.theme() +
facet_wrap(~site) +
theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
strip.background = element_rect(fill = "lightgrey"),
axis.title = element_text(size = 26),
axis.text = element_text(size = 26),
legend.text = element_text(size = 26),
legend.position = "bottom",
strip.text.x = element_text(size = 26)) +
scale_y_continuous(limits = c(0,25), breaks = seq(0, 25, by = 5),
name = "Species Richness") +
scale_x_continuous(limits = c(0,28), breaks = seq(0, 30, by = 7),
name = "Survey Effort (Days)") +
scale_colour_manual(values = c("#FF8C00", "#A034F0", "#189F9F","#D14103","black"),
name = "",
labels = c("OBM","Camera Trap", "PAM (survey period)", "Remote Sensing",
"Total Richness")))
Vocal Mammals
Considering vocal mammals only and excluding non-vocal mammals from the analysis.
Prepare data
Merged duplicate rows for duplicate observations and transformed the dataset to a wide format.
mammal_community_vocal <- mammal_data_vocal %>%
select(site, SamplingDay, assessment.method,common.name) %>% #selected relevant grouping variables
group_by(across(everything())) %>% #groups by all available columns
mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
ungroup() %>%
unique() %>% #merged the duplicate rows
pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
values_from = "number",
values_fill = list(number=0)) %>% #replaced NAs with 0
arrange(site) %>%
filter(assessment.method != "PAM.all") %>%
mutate(SamplingDay = as.factor(SamplingDay))Added a total category to assessment methods.
result.i <- vector("list",length(unique(mammal_community_vocal$site))) #created an empty list that was filled as the loop ran
for(i in 1:length(unique(mammal_community_vocal$site))){ #looped over each study site
site.i <- mammal_community_vocal[mammal_community_vocal$site == unique(mammal_community_vocal$site)[i],] #subset to the site
result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] #subset to the SamplingDay
meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("total"))
SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
colnames(total.a) <- colnames(mammal_community_vocal)
result.a[[a]] <- total.a} #added method results to list
result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop
result.i <- do.call("rbind.data.frame",result.i)
mammal_community_vocal2 <- rbind.data.frame(mammal_community_vocal, result.i)Added a remote sensing category to assessment methods. Remote sensing was a combination of PAM and camera trapping.
result.j <- vector("list",length(unique(mammal_community_vocal$site))) #created an empty list that was filled as the loop ran
for(i in 1:length(unique(mammal_community_vocal$site))){ #looped over each study site
site.i <- mammal_community_vocal[mammal_community_vocal$site == unique(mammal_community_vocal$site)[i],] #subset to the site
result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] %>%
filter(assessment.method %in% c("camera", "PAM.survey"))
meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("remote"))
SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
colnames(total.a) <- colnames(mammal_community_vocal)
result.a[[a]] <- total.a} #added method results to list
result.j[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop
result.j <- do.call("rbind.data.frame",result.j)
mammal_community_vocal2 <- rbind.data.frame(mammal_community_vocal2, result.j) %>%
drop_na()Loop
Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.
# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 30, ncol = (ncol(mammal_community_vocal2)))) %>%
replace(is.na(.), 0) %>%
set_names(names(mammal_community_vocal2)) %>%
select(-(1:3))
result.i <- vector("list",length(unique(mammal_community_vocal2$site))) #created an empty list that was filled as the loop ran
set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(mammal_community_vocal2$site))){ #looped over each study site
site.i <- mammal_community_vocal2[mammal_community_vocal2$site == unique(mammal_community_vocal2$site)[i],] #subset to the site
dates.i <- ifelse(unique(site.i$site) %in% c("Tarcutta", "Undara", "Wambiana", "Rinyirru"), 28, 21) #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day
result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
sampling.a <- method.a$assessment.method[1] #saved the name of the method
method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
accum.a <- specaccum(method.a,"random") #calculated data
#created data frame of metadata, richness, and standard deviation
accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
result.a[[a]] <- accum.a}} #added method results to list
result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loopTransformed resulting list to a data frame.
mammals_wide.result_vocal <- do.call("rbind.data.frame",result.i) %>% #compressed result.i list into data frame
mutate(assessment.method = factor(assessment.method, levels = c("OBM",
"camera",
"PAM.survey",
"remote",
"total")),
site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))Visualisation
Plot of rarefaction curves for each method split by sites.
(specaccum_mammals_vocal <- ggplot(mammals_wide.result_vocal, aes(x = day, y = richness, col = assessment.method)) +
geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.3))),
linetype = 0.1) +
geom_vline(xintercept=c(7,14,21, 28), linetype="dotted", colour = "black") +
geom_line(aes(group = assessment.method), linewidth = 1.3) +
my.theme() +
facet_wrap(~site) +
theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
strip.background = element_rect(fill = "lightgrey"),
axis.title = element_text(size = 26),
axis.text = element_text(size = 26),
legend.text = element_text(size = 26),
legend.position = "bottom",
strip.text.x = element_text(size = 26)) +
scale_y_continuous(limits = c(0,15), breaks = seq(0, 15, by = 5),
name = "") +
scale_x_continuous(limits = c(0,28), breaks = seq(0, 30, by = 7),
name = "Survey Effort (Days)") +
scale_colour_manual(values = c("#FF8C00", "#A034F0", "#189F9F","#D14103" ,"black"),
name = "",
labels = c("OBM","Camera Trap", "PAM (survey period)", "Remote Sensing",
"Total Richness")))
Combined Figure
(specaccum_mammals_combined <- (specaccum_mammals + ggtitle('All mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(specaccum_mammals_vocal + ggtitle('Vocal mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
theme(legend.position = c(-0.1,-0.2), legend.text = element_text(size = 22),
plot.margin = unit(c(5, 10, 30, 10), units = "mm")) +
guides(colour = guide_legend(nrow = 1)))
PAM survey length
How many days until PAM.all reaches values similar to OBM for all mammals and only vocal mammals. In other words how long should we leave recorders out on average until we get the best number of species.
all_mammals <- mammal_data %>%
filter(assessment.method == "PAM.all") %>%
select(site, SamplingDay, assessment.method,common.name) %>% #selected relevant grouping variables
group_by(across(everything())) %>% #groups by all available columns
mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
ungroup() %>%
unique() %>% #merged the duplicate rows
pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
values_from = "number",
values_fill = list(number=0)) %>% #replaced NAs with 0
arrange(site) %>%
mutate(SamplingDay = as.factor(SamplingDay))Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.
# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 1235, ncol = (ncol(all_mammals)))) %>%
replace(is.na(.), 0) %>%
set_names(names(all_mammals)) %>%
select(-(1:3))
result.i <- vector("list",length(unique(all_mammals$site))) #created an empty list that was filled as the loop ran
set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(all_mammals$site))){ #looped over each study site
site.i <- all_mammals[all_mammals$site == unique(all_mammals$site)[i],] #subset to the site
#manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day
site_values <- c(Tarcutta = 986, Duval = 679, Mourachan = 783, Wambiana = 1235, Undara = 780, Rinyirru = 488)
dates.i <- site_values[match(unique(site.i$site), names(site_values))]
result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
sampling.a <- method.a$assessment.method[1] #saved the name of the method
method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
accum.a <- specaccum(method.a,"random") #calculated data
#created data frame of metadata, richness, and standard deviation
accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
result.a[[a]] <- accum.a}} #added method results to list
result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loopTransformed resulting list to a data frame.
all_mammals.result <- do.call("rbind.data.frame",result.i) %>% #compressed result.i list into data
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))pam.max <- all_mammals.result %>%
group_by(site) %>%
filter(assessment.method == "PAM.all",
day == max(day)) %>%
ungroup() %>%
select(site, richness.max = richness)
pam.28 <- all_mammals.result %>%
group_by(site) %>%
filter(assessment.method == "PAM.all",
day == 28) %>%
ungroup() %>%
select(richness.28 = richness) %>%
mutate(prop.28 = round(richness.28/pam.max$richness.max, digits = 2))
pam.90 <- all_mammals.result %>%
group_by(site) %>%
filter(assessment.method == "PAM.all",
day == 90) %>%
ungroup() %>%
select(richness.90 = richness) %>%
mutate(prop.90 = round(richness.90/pam.max$richness.max, digits = 2))
pam.180 <- all_mammals.result %>%
group_by(site) %>%
filter(assessment.method == "PAM.all",
day == 180) %>%
ungroup() %>%
select(richness.180 = richness) %>%
mutate(prop.180 = round(richness.180/pam.max$richness.max, digits = 2))
pam.365 <- all_mammals.result %>%
group_by(site) %>%
filter(assessment.method == "PAM.all",
day == 365) %>%
ungroup() %>%
select(richness.365 = richness) %>%
mutate(prop.365 = round(richness.365/pam.max$richness.max, digits = 2))
pam.rich.prop <- cbind(pam.max,pam.28,pam.90,pam.180,pam.365) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
(sampling.table <- pam.rich.prop %>%
select(Site = site, "PAM (total richness)" = richness.max,
"PAM (proportion 28 days)" = prop.28,
"PAM (proportion 90 days)" = prop.90,
"PAM (proportion 180 days)" = prop.180,
"PAM (proportion 365 days)" = prop.365) %>%
add_row(Site = "Average", "PAM (total richness)" = round(mean(pam.rich.prop$richness.max)),
"PAM (proportion 28 days)" = round(mean(pam.rich.prop$prop.28), digits = 2),
"PAM (proportion 90 days)" = round(mean(pam.rich.prop$prop.90), digits = 2),
"PAM (proportion 180 days)" = round(mean(pam.rich.prop$prop.180), digits = 2),
"PAM (proportion 365 days)" = round(mean(pam.rich.prop$prop.365), digits = 2))) Site PAM (total richness) PAM (proportion 28 days) PAM (proportion 90 days) PAM (proportion 180 days)
1 Duval 10 0.72 0.86 0.91
2 Mourachan 7 0.68 0.86 0.89
3 Rinyirru 9 0.61 0.79 0.87
4 Tarcutta 13 0.66 0.78 0.86
5 Undara 7 0.50 0.75 0.86
6 Wambiana 10 0.79 0.91 0.96
7 Average 9 0.66 0.82 0.89
PAM (proportion 365 days)
1 0.95
2 0.92
3 0.96
4 0.94
5 0.92
6 0.99
7 0.95Visualisation
vline_data <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
xintercept = c(12, 175, 17, 88, 676,320)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
vline_data2 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
xintercept = c(12, 35, 12, 17, 12,9)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
hline_data_vocal <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
yintercept = c(4.1, 6, 7, 6,12,10)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
hline_data_vocal2 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
yintercept = c(3.9, 4, 6, 4,6,6)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
hline_data_cam <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
yintercept = c(4, 4, 6, 4,6,6)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
vline_data_cam <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
xintercept = c(12, 35, 12, 17, 12,9)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
hline_data_28 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
yintercept = c(5.48, 3.51, 7.88, 4.32,8.52,6.87)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
vline_data_28 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
xintercept = c(28, 28, 28, 21, 21,28)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
hline_data_total <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
yintercept = c(8, 7, 11, 8,13,11)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))
vline_data_total <- data.frame(site = c("Rinyirru", "Undara","Tarcutta"),
xintercept = c(210, 771, 968)) %>%
mutate(site = factor(site, levels = c("Rinyirru", "Undara","Tarcutta")))
(spec_all_mammals_plot <- ggplot(all_mammals.result, aes(x = day, y = richness, col = assessment.method)) +
geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.5))),
linetype = 0.1) +
geom_line(aes(group = assessment.method), linewidth = 1.5) +
my.theme() +
facet_wrap(~site, scales = "free_x") +
geom_hline(data = hline_data_vocal, aes(yintercept = yintercept), linetype="dashed", color = "#FF8C00",
linewidth = 1) +
geom_hline(data = hline_data_vocal2, aes(yintercept = yintercept), linetype="dashed", color = "#A034F0",
linewidth = 1) +
#geom_vline(data = vline_data, aes(xintercept = xintercept), linetype="dotted", color = "#FF8C00") +
#geom_vline(data = vline_data2, aes(xintercept = xintercept), linetype="dotted", color = "#A034F0") +
geom_text(data = vline_data, aes(x = 100, y = 14, label = paste0("DTM: ",xintercept)),
vjust = -1, size = 8,inherit.aes = FALSE, color = "#FF8C00") +
geom_text(data = vline_data2, aes(x = 100, y = 13, label = paste0("DTM: ",xintercept)),
vjust = -1, size = 8,inherit.aes = FALSE, color = "#A034F0") +
theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
strip.background = element_rect(fill = "lightgrey"),
axis.title = element_text(size = 26),
axis.text = element_text(size = 26),
legend.text = element_text(size = 26),
legend.position = "bottom",
strip.text.x = element_text(size = 26)) +
scale_y_continuous(limits = c(0,15), breaks = seq(0, 15, by = 5),
name = "Species Richness") +
scale_x_continuous(name = "Survey Effort (Days)") +
scale_colour_manual(values = c("#189F9F"),
name = "",
labels = c("PAM (all audio)")))
Session Info
sessionInfo()R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Australia/Brisbane
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] pilot_4.0.1 paletteer_1.6.0 cowplot_1.1.3 ggimage_0.3.3 ggsignif_0.6.4
[6] report_0.6.1 HDInterval_0.2.4 patchwork_1.3.0 EcolUtils_0.1 vegan_2.6-10
[11] lattice_0.22-6 permute_0.9-7 tidybayes_3.0.7 emmeans_1.11.1 DHARMa_0.4.7
[16] ggeffects_2.2.1 rstan_2.32.7 StanHeaders_2.32.10 brms_2.22.0 Rcpp_1.0.14
[21] cmdstanr_0.9.0 colorspace_2.1-1 lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1
[26] dplyr_1.1.4 purrr_1.0.4 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[31] ggplot2_3.5.2 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 tensorA_0.36.2.1 rstudioapi_0.17.1 jsonlite_2.0.0 magrittr_2.0.3
[6] estimability_1.5.1 magick_2.8.6 farver_2.1.2 nloptr_2.2.1 rmarkdown_2.29
[11] fs_1.6.6 vctrs_0.6.5 minqa_1.2.8 htmltools_0.5.8.1 distributional_0.5.0
[16] curl_6.2.2 gridGraphics_0.5-1 htmlwidgets_1.6.4 plyr_1.8.9 mime_0.13
[21] iterators_1.0.14 lifecycle_1.0.4 pkgconfig_2.0.3 gap_1.6 Matrix_1.7-3
[26] R6_2.6.1 fastmap_1.2.0 shiny_1.10.0 rbibutils_2.3 digest_0.6.37
[31] rematch2_2.1.2 ps_1.9.1 qgam_2.0.0 labeling_0.4.3 timechange_0.3.0
[36] abind_1.4-8 mgcv_1.9-1 compiler_4.5.0 doParallel_1.0.17 bit64_4.6.0-1
[41] withr_3.0.2 backports_1.5.0 inline_0.3.21 QuickJSR_1.7.0 pkgbuild_1.4.7
[46] MASS_7.3-65 loo_2.8.0 tools_4.5.0 httpuv_1.6.16 glue_1.8.0
[51] promises_1.3.2 nlme_3.1-168 grid_4.5.0 checkmate_2.3.2 cluster_2.1.8.1
[56] reshape2_1.4.4 generics_0.1.4 gtable_0.3.6 tzdb_0.5.0 data.table_1.17.2
[61] hms_1.1.3 utf8_1.2.5 foreach_1.5.2 pillar_1.10.2 ggdist_3.3.3
[66] yulab.utils_0.2.0 vroom_1.6.5 later_1.4.2 posterior_1.6.1 splines_4.5.0
[71] bit_4.6.0 tidyselect_1.2.1 knitr_1.50 reformulas_0.4.1 arrayhelpers_1.1-0
[76] gridExtra_2.3 V8_6.0.5 stats4_4.5.0 xfun_0.52 bridgesampling_1.1-2
[81] matrixStats_1.5.0 stringi_1.8.7 ggfun_0.1.8 yaml_2.3.10 boot_1.3-31
[86] evaluate_1.0.3 codetools_0.2-20 ggplotify_0.1.2 cli_3.6.5 RcppParallel_5.1.10
[91] xtable_1.8-4 Rdpack_2.6.4 processx_3.8.6 coda_0.19-4.1 svUnit_1.0.6
[96] parallel_4.5.0 rstantools_2.4.0 gap.datasets_0.0.6 bayesplot_1.12.0 Brobdingnag_1.2-9
[101] lme4_1.1-37 mvtnorm_1.3-3 scales_1.4.0 insight_1.3.1 crayon_1.5.3
[106] rlang_1.1.6 Cite Packages
cite_packages() - Bürkner P (2017). "brms: An R Package for Bayesian Multilevel Models Using Stan." _Journal of Statistical Software_, *80*(1), 1-28. doi:10.18637/jss.v080.i01 <https://doi.org/10.18637/jss.v080.i01>. Bürkner P (2018). "Advanced Bayesian Multilevel Modeling with the R Package brms." _The R Journal_, *10*(1), 395-411. doi:10.32614/RJ-2018-017 <https://doi.org/10.32614/RJ-2018-017>. Bürkner P (2021). "Bayesian Item Response Modeling in R with brms and Stan." _Journal of Statistical Software_, *100*(5), 1-54. doi:10.18637/jss.v100.i05 <https://doi.org/10.18637/jss.v100.i05>.
- Constantin A, Patil I (2021). "ggsignif: R Package for Displaying Significance Brackets for 'ggplot2'." _PsyArxiv_. doi:10.31234/osf.io/7awm6 <https://doi.org/10.31234/osf.io/7awm6>, <https://psyarxiv.com/7awm6>.
- Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2025). _Rcpp: Seamless R and C++ Integration_. doi:10.32614/CRAN.package.Rcpp <https://doi.org/10.32614/CRAN.package.Rcpp>, R package version 1.0.14, <https://CRAN.R-project.org/package=Rcpp>. Eddelbuettel D, François R (2011). "Rcpp: Seamless R and C++ Integration." _Journal of Statistical Software_, *40*(8), 1-18. doi:10.18637/jss.v040.i08 <https://doi.org/10.18637/jss.v040.i08>. Eddelbuettel D (2013). _Seamless R and C++ Integration with Rcpp_. Springer, New York. doi:10.1007/978-1-4614-6868-4 <https://doi.org/10.1007/978-1-4614-6868-4>, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). "Extending R with C++: A Brief Introduction to Rcpp." _The American Statistician_, *72*(1), 28-36. doi:10.1080/00031305.2017.1375990 <https://doi.org/10.1080/00031305.2017.1375990>.
- Gabry J, Češnovar R, Johnson A, Bronder S (2025). _cmdstanr: R Interface to 'CmdStan'_. R package version 0.9.0, <https://mc-stan.org/cmdstanr/>.
- Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate." _Journal of Statistical Software_, *40*(3), 1-25. <https://www.jstatsoft.org/v40/i03/>.
- Hartig F (2024). _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models_. doi:10.32614/CRAN.package.DHARMa <https://doi.org/10.32614/CRAN.package.DHARMa>, R package version 0.4.7, <https://CRAN.R-project.org/package=DHARMa>.
- Hawkins O (2025). _pilot: A minimal ggplot2 theme with an accessible discrete color palette_. R package version 4.0.1, commit 532e68c5135fc847130a76547f3d32832bf43f08, <https://github.com/olihawkins/pilot>.
- Hvitfeldt E (2021). _paletteer: Comprehensive Collection of Color Palettes_. R package version 1.3.0, <https://github.com/EmilHvitfeldt/paletteer>.
- Kay M (2024). _tidybayes: Tidy Data and Geoms for Bayesian Models_. doi:10.5281/zenodo.1308151 <https://doi.org/10.5281/zenodo.1308151>, R package version 3.0.7, <http://mjskay.github.io/tidybayes/>.
- Lenth R (2025). _emmeans: Estimated Marginal Means, aka Least-Squares Means_. doi:10.32614/CRAN.package.emmeans <https://doi.org/10.32614/CRAN.package.emmeans>, R package version 1.11.1, <https://CRAN.R-project.org/package=emmeans>.
- Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects from Regression Models." _Journal of Open Source Software_, *3*(26), 772. doi:10.21105/joss.00772 <https://doi.org/10.21105/joss.00772>.
- Makowski D, Lüdecke D, Patil I, Thériault R, Ben-Shachar M, Wiernik B (2023). "Automated Results Reporting as a Practical Tool to Improve Reproducibility and Methodological Best Practices Adoption." _CRAN_. <https://easystats.github.io/report/>.
- Meredith M, Kruschke J (2022). _HDInterval: Highest (Posterior) Density Intervals_. doi:10.32614/CRAN.package.HDInterval <https://doi.org/10.32614/CRAN.package.HDInterval>, R package version 0.2.4, <https://CRAN.R-project.org/package=HDInterval>.
- Müller K, Wickham H (2023). _tibble: Simple Data Frames_. doi:10.32614/CRAN.package.tibble <https://doi.org/10.32614/CRAN.package.tibble>, R package version 3.2.1, <https://CRAN.R-project.org/package=tibble>.
- Oksanen J, Simpson G, Blanchet F, Kindt R, Legendre P, Minchin P, O'Hara R, Solymos P, Stevens M, Szoecs E, Wagner H, Barbour M, Bedward M, Bolker B, Borcard D, Carvalho G, Chirico M, De Caceres M, Durand S, Evangelista H, FitzJohn R, Friendly M, Furneaux B, Hannigan G, Hill M, Lahti L, McGlinn D, Ouellette M, Ribeiro Cunha E, Smith T, Stier A, Ter Braak C, Weedon J, Borman T (2025). _vegan: Community Ecology Package_. doi:10.32614/CRAN.package.vegan <https://doi.org/10.32614/CRAN.package.vegan>, R package version 2.6-10, <https://CRAN.R-project.org/package=vegan>.
- Pedersen T (2024). _patchwork: The Composer of Plots_. doi:10.32614/CRAN.package.patchwork <https://doi.org/10.32614/CRAN.package.patchwork>, R package version 1.3.0, <https://CRAN.R-project.org/package=patchwork>.
- R Core Team (2025). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.
- Salazar G (2025). _EcolUtils: Utilities for community ecology analysis_. R package version 0.1, commit 2b57a8a8b76eba0852c56ce30f93d92c30320161, <https://github.com/GuillemSalazar/EcolUtils>.
- Sarkar D (2008). _Lattice: Multivariate Data Visualization with R_. Springer, New York. ISBN 978-0-387-75968-5, <http://lmdvr.r-forge.r-project.org>.
- Simpson G (2022). _permute: Functions for Generating Restricted Permutations of Data_. doi:10.32614/CRAN.package.permute <https://doi.org/10.32614/CRAN.package.permute>, R package version 0.9-7, <https://CRAN.R-project.org/package=permute>.
- Stan Development Team (2020). "StanHeaders: Headers for the R interface to Stan." R package version 2.21.0-6, <https://mc-stan.org/>.
- Stan Development Team (2025). "RStan: the R interface to Stan." R package version 2.32.7, <https://mc-stan.org/>.
- Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_. Springer-Verlag New York. ISBN 978-3-319-24277-4, <https://ggplot2.tidyverse.org>.
- Wickham H (2023). _forcats: Tools for Working with Categorical Variables (Factors)_. doi:10.32614/CRAN.package.forcats <https://doi.org/10.32614/CRAN.package.forcats>, R package version 1.0.0, <https://CRAN.R-project.org/package=forcats>.
- Wickham H (2023). _stringr: Simple, Consistent Wrappers for Common String Operations_. doi:10.32614/CRAN.package.stringr <https://doi.org/10.32614/CRAN.package.stringr>, R package version 1.5.1, <https://CRAN.R-project.org/package=stringr>.
- Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." _Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
- Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar of Data Manipulation_. doi:10.32614/CRAN.package.dplyr <https://doi.org/10.32614/CRAN.package.dplyr>, R package version 1.1.4, <https://CRAN.R-project.org/package=dplyr>.
- Wickham H, Henry L (2025). _purrr: Functional Programming Tools_. doi:10.32614/CRAN.package.purrr <https://doi.org/10.32614/CRAN.package.purrr>, R package version 1.0.4, <https://CRAN.R-project.org/package=purrr>.
- Wickham H, Hester J, Bryan J (2024). _readr: Read Rectangular Text Data_. doi:10.32614/CRAN.package.readr <https://doi.org/10.32614/CRAN.package.readr>, R package version 2.1.5, <https://CRAN.R-project.org/package=readr>.
- Wickham H, Vaughan D, Girlich M (2024). _tidyr: Tidy Messy Data_. doi:10.32614/CRAN.package.tidyr <https://doi.org/10.32614/CRAN.package.tidyr>, R package version 1.3.1, <https://CRAN.R-project.org/package=tidyr>.
- Wilke C (2024). _cowplot: Streamlined Plot Theme and Plot Annotations for 'ggplot2'_. doi:10.32614/CRAN.package.cowplot <https://doi.org/10.32614/CRAN.package.cowplot>, R package version 1.1.3, <https://CRAN.R-project.org/package=cowplot>.
- Yu G (2023). _ggimage: Use Image in 'ggplot2'_. doi:10.32614/CRAN.package.ggimage <https://doi.org/10.32614/CRAN.package.ggimage>, R package version 0.3.3, <https://CRAN.R-project.org/package=ggimage>.
- Zeileis A, Fisher JC, Hornik K, Ihaka R, McWhite CD, Murrell P, Stauffer R, Wilke CO (2020). "colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes." _Journal of Statistical Software_, *96*(1), 1-49. doi:10.18637/jss.v096.i01 <https://doi.org/10.18637/jss.v096.i01>. Zeileis A, Hornik K, Murrell P (2009). "Escaping RGBland: Selecting Colors for Statistical Graphics." _Computational Statistics & Data Analysis_, *53*(9), 3259-3270. doi:10.1016/j.csda.2008.11.033 <https://doi.org/10.1016/j.csda.2008.11.033>. Stauffer R, Mayr GJ, Dabernig M, Zeileis A (2009). "Somewhere over the Rainbow: How to Make Effective Use of Colors in Meteorological Visualizations." _Bulletin of the American Meteorological Society_, *96*(2), 203-216. doi:10.1175/BAMS-D-13-00155.1 <https://doi.org/10.1175/BAMS-D-13-00155.1>.